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Unintegrated gluon distributions sensitive to the transverse spatial distribution of gluons in the 
proton are extracted from data on exclusive and diffractive final states at HERA in the dipole 
approach. These unintegrated gluon distributions can be used to compute inclusive hadron pro- 
duction in p+p collisions at the LHC. In this paper, we consider a number of saturation models 
with differing dynamical assumptions that give good fits to the available HERA data. We apply 
these models to study the rapidity and transverse momentum dependence of the LHC data up to 
*Js — 7 TeV. We examine the sensitivity of these results to parameters that are not constrained by 
the HERA data and comment on similarities and differences with previous work. We compute the 
n-particle inclusive multiplicity distribution and show that the LHC p+p results are in agreement 
with predictions for multi-particle production in the Color Glass Condensate approach. This result 
has significant ramifications for the interpretation of multi-particle correlations in high multiplicity 
events at the LHC. 

I. INTRODUCTION 

Colliders at ultra-high energies have made feasible the study of the collective many-body properties of QCD in the 
limit of fixed Q 2 ^ Aq CD and very small x. An important property of QCD in these "Regge asymptotics" is the 
phenomenon of gluon saturation, which states that a probe with transverse resolution 1/Q 2 interacts with a target 
with probability of order unity when Q 2 < Q 2 s (x), where Q 2 s {x) is a dynamically generated semi- hard scale in the 
target hadron or nucleus In the framework where the parton model is manifest, the gluon saturation phenomenon 
corresponds to the requirement that transverse momentum modes with k± < Q$ are maximally occupied with an 
occupation number 1/as- The large occupation number means that these modes can be treated classically Q; an 
additional separation of scale between fast and slow modes due to time dilation allows for a renormalization group 
treatment of these high occupation number states as they evolve with energy Q. This Color Glass Condensate 
(CGC) [4| description of saturated gluons is universal to hadrons and nuclei and can be tested both in deeply inelastic 
00 , scattering experiments off hadrons and nuclei as well as in hadron/nuclear collisions. 

The CGC approach, when applied to Deeply Inelastic Scattering (DIS), results at leading order in a s [gj in the 
dipole picture where the inclusive virtual photon hadron cross section is expressed as @ 
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Here 
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*X r (r ± ,z,Q) 



represents the probability for a virtual photon to produce a quark-anti-quark pair of size 

r = |rj_| and d i^' p (t±, x, h±) denotes the dipole cross section for this pair to scatter off the target at an impact 
parameter b^. Tne former is well known from QED, while the latter represents the dynamics of QCD scattering 
at small x. It was shown some time ago that a model known as the Golec-Biernat-Wusthoff (GBW) model [7] that 

implements saturation in the dipole cross-section through the parametrization = 2(1 — e~ r Qs.p^)/ 4 ), where 

Qli-p{ x ) — ( x o/ x ) X GeV 2 , gives a good qualitative fit to the HERA inclusive and diffractive cross section data for 
x n = 3 • 10 -4 and A = 0.288. This work explained very simply key features of the HERA data and was very suggestive 
of the possible role of a semi-hard saturation scale in the hadron. This work was further refined in more sophisticated 
models that are somewhat better motivated and treat the impact parameter dependence of the dipole cross-section 
more accurately. As we shall discuss further in the next section, these models give excellent fits to small x inclusive, 
diffractive and exclusive HERA data. The common ingredient in these combined fits is the dipole cross-section. 

The dipole cross-section, to leading logarithmic accuracy, is a universal quantity which can be applied to compute 
inclusive quantities in hadron-hadron collisions. It is defined in terms of the real part of the forward scattering 
amplitude Af(r±,x, bj_) as 

Un " ip -(r x ,x,b ± ) = 2A^(r ± ,x,b ± )^2 fl-^(tr(^(b ± + ^)?7t ( b ± -^))\ ) , (2) 



2 





FIG. 1: Left: Dipole cross-section in DIS. Right: Overlap of unintegrated gluon distributions in proton-proton collisions. 



where U(h± ± is a Wilson line in the fundamental representation representing the interaction between a quark 
and the color fields of the target. The average (• • • ) x is an average over these color fields; the energy dependence of 
the correlator as a function of x (or the rapidity Y — \sx{l/x)) is given by the JIMWLK equation j3j. In the large N c 
limit, the equation for the energy evolution of this correlator is the Balitsky-Kovchegov (BK) equation Q. We note 
however that neither JIMWLK nor BK is at present equipped to deal well with the impact parameter dependence of 
the dipole cross-section; the dipole cross-section in this formalism is taken in eq. ^ to be independent of the impact 
parameter. To address the impact parameter dependence of this equation, one resorts to models which parametrize 
both saturation effects and the impact parameter dependence. 

In hadron-hadron collisions, one can derive at leading order the expression [9| 



dN g (h ± ) I6a s 1 



d 2 k 



dV 



d(/>A(xi,kj_|s_L) d0 B (>2,P± - kjs ± - bj 



This equation is a generalization of the well known k±_ factorization expression for inclusive gluon production [Tpj 
to include the impact parameter dependence of the unintegrated gluon distributions. Here Cf = — 1/2N C is the 
Casimir for the fundamental representation. Using a relation between quark and gluon dipole amplitudes strictly 
valid in the large N c limit, the unintegrated gluon distribution in either of the two protons can be expressed in terms 
of the corresponding dipole cross-section measured in DIS as (Tlj 
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(4) 



Thus the impact parameter dependent dipole cross-section determined from HERA data can be used to compute 
the single inclusive gluon distribution in proton-proton collisions with no additional parameters. This statement is 
strictly valid to leading log accuracy for momenta k± > Q 2 p . Further, as we shall discuss later, there will be additional 
parameters that come in when one wants to make contact with the measured hadron spectrum. 

This approach was applied most recently to compute the single inclusive hadron spectrum in proton-proton collisions 
at the LHC by Levin and Rezaeian pjj. The quantitative differences of our study to their work are the following: 
a) we consider three dipole models that give good fits to HERA data to see whether they give results consistent 
with the LHC data, b) we study and comment on the dependence of the results on variations of the parameters in 
the study and c) we convolve the inclusive gluon distribution with fragmentation function instead of using a simple 
fragmentation presciption as in ref. [12J. We shall also comment on other quantitative differences in our respective 
treatments at the appropriate points in the text. A qualitative difference of our work relative to that of ref. [l2| is that 
we compute directly the average inclusive multiplicity at a given impact parameter. In computing the minimum bias 
single inclusive multiplicity distribution, there are similar uncertainties as ref. [l2| , which can be fixed by normalizing 
the data to single inclusive data at lower center of mass energies. However, as we shall discuss later, the average 
multiplicity at a given impact parameter is an essential input in computing the probability distribution as a function 
of event multiplicity. We shall compute the n-particle probability distribution and compare our results with the p+p 
collider data. These results will be important in understanding the role of various sources of fluctuations in the p+p 
collider data. 
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II. PROTON DIPOLE CROSS-SECTION 



In this section, we shall discuss some of the dipole model parametrizations that have been compared extensively to 
the HERA data. We shall consider respectively the IP-Sat model [16], the b-CGC model [l3l4Ta | and the NLO-BK 
model [13] • Our list of dipole parametrizations is by no means exhaustive [l8[ but is a sample of some of the key 
approaches where extensive comparisons have been made to the data. We note that dipole models have been applied 
to understand RHIC data Il9f . and there have also been preliminary attempts at a combined analysis of RHIC, HERA 
and fixed target e+ A data [20-24] . We also note that there are leading twist models with some similarities to dipole 
models that have been compared to the HERA data (25j . 



The IP-Sat Model 



The impact parameter dependent dipole saturation model (IP-Sat) is a refinement of the Golec-Biernat- 
Wusthoff dipole model Q to give the right perturbative limit when r_L — > [2f|. It is equivalent to the expression 
derived in the classical effective theory of the CGC, to leading logarithmic accuracy @. The proton dipole cross-section 
in this model is expressed as 



da d p ip 



(rx.ar.bj.) = 2 



d 2 b ± 

Here the scale /x 2 is related to dipole radius r± (see fig. [J) as 

4 
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where the leading order expression for the running coupling is 



as(/J 2 ) = 



12tt 



(33-2n / )log(/ J 2 /A 2 3CD ) 



(5) 



(6) 



(7) 



with n/ =3, Aqcd=0.2 GeV. The model includes saturation as eikonalized power corrections to the DGLAP leading 
twist expression and may be valid in the regime where logs in Q 2 dominate logs in x. The saturation scale for a 
fixed impact parameter is determined self-consistently by requiring that the dipole amplitude (within brackets in 
cq. [5]) have the magnitude Af(x, rs, bj_) = 1 — e -1 / 2 , with Q 2 p — 2/r|. We note that there is an overall logarithmic 
uncertainty in the determination of Q 2 p (x, b±). 

For each value of the dipole radius, the gluon density xg(x, fj, 2 ) is evolved from //§ to /i 2 using LO DGLAP evolution 
equation without quarks, 



dxg(x,n 2 ) _ as{n 2 ) 
d log /j 2 2tt 



X fX n\ 

dzP gg( z )-9 ) 



(8) 



Here the gluon splitting function with 71/ flavor and CU=3 & Tr—1 is 

z 1 — z 



PggW = 6 



+ z(l-z) 



.(!"*)+ 

The initial gluon density at the scale (1q is taken to be of the form 

xg(x, /Jq) = A g x~ x <> (1 - a;) 5 - 6 



(9) 



(10) 



An important feature of the IP-Sat model is the b-dependence of the dipole cross-section, which is introduced through 
a gluon density profile function T(b). This profile function is normalized to unity and is chosen to have the Gaussian 
form 



T P (b ± ) = 



1 



2nB G 



exp 



-bL a 
2B G 



where Bq is a parameter fit to the HERA diffractive data. This corresponds to (b 2 ) 
gluonic radius of the proton. 



(11) 

2Bq, the average squared 
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B G (GeV" 2 ) 


MGeV 2 ) 


A a 


A 9 


1.4 


4.0 


1.17 


2.55 


0.020 


1.35 


4.0 


1.20 


2.51 


0.024 


1.5 


4.0 


0.77 


2.64 


0.011 


1.4 


4.0 


1.50 


3.61 


-0.118 



TABLE I: Parameters of the IP-Sat model obtained from the fit to HERA data Q. 



Sets of parameters obtained from optimal fits of the IP-Sat model to HERA data 14j are listed in table HI All data 
sets except the last use m Uj d, s = 0.14 GeV; the last set corresponds to m„ a a = 0.05 GeV. The parameters of the 
initial gluon distribution are determined from fits to the HERA data [23, [28( with a \ 2 ~ 1- For charm quarks, 
x = Xbj(l + 4m 2 /Q 2 ). The value of Bq is determined primarily from the ^-distributions of J/ip mesons measured 
by ZEUS [29J and HI [30j. With these parameters, excellent agreement is obtained with the HERA exclusive vector 
meson and DVCS data. For a detailed comparison of this model to the HERA data, we refer the reader to Ref. [lH • 

B. The b-CGC Model 

The IP-sat dipole model is applicable when leading logarithms in Q 2 dominate over leading logarithms in x. At 
very small x, quantum evolution in the CGC describing both the bremsstrahlung limit of linear small x evolution 
as well as nonlinear RG evolution at high parton densities, combined with a realistic 6-dependence, may be better 
captured in the bCGC model [p^ - fTo] ]. The proton dipole cross-section in this case is expressed as 

- 5 £P(r L ,x J b_ L ) =2x{ ■»<> {—) ■ r ^ ^ 2 > (12) 

^ [ 1 - exp (-A ln 2 (Srj.Q s )) : r ± Q s > 2; 

In this model, in contrast to the IP-sat model, the impact parameter dependence is introduced though the quantity 
Q s (x,b±), denned as 



exp 



2B CGC ) 



(13) 



As previously, for comparison of scales among different saturation models, the relevant saturation scale for a fixed 
impact parameter is determined self-consistently by requiring that the dipole amplitude (the expression to the right 
of the curly bracket in eq. IT2"j) have the magnitude Af(x, rs, bjj = 1 — e -1 / 2 , with the saturation scale denned as 
Q 2 p = 2/r|. The coefficients A and B are obtained by requiring the two asymptotic forms of the dipole cross-section 
and their first derivatives are continuous at r±Q s = 2: 

A = d-w-Mi k B = 5 (1 - «r^ («) 

The parameter k = 9.9 is fixed from the leading order BFKL value for this quantity 

Table |H] presents the parameters of the model that are fitted to the HERA data . The parameter Bqgc is 
determined from the t-dependence of exclusive J/tp photo-production. For the b-CGC model, this parameter however 
cannot be easily interpreted as giving the square mean gluonic radius of the proton. The parameters presented in the 
table in the first and fourth lines do not give good fits to the data. The fit corresponding to the second line of the table 
gives the best fit to the data with % 2 ~ 1. The third line of the table corresponds to a fit where no saturation form is 
employed (namely, the perturbative expression, without the diffusion term proportional to Y = ln(xo/x), is extended 
to r±Q s > 2); it gives equally good fits to the data. However, it should be noted that this choice of parameters will 
violate perturbative unitarity for large dipole sizes > 1/Q S . 

The dependence of the saturation scale Q 2 p (x,bj^) as a function of x for different (and vice versa) in the 
IP-Sat and b-CGC models is shown in fig. [5J In both cases, the fits to the HERA data result in a semi-hard scale 
{Qs ^QCd) w itli decreasing x and b values probed in the collisions. The existence of such scales and their increase 
with energy is what validates the whole approach of treating high parton densities in weak coupling. It would of 
course be naive to interpret the extracted numerical value of Qs as being precisely the scale that controls the running 
of the coupling. As is well known, the scale that controls the running of the coupling can differ considerably from this 
"bare" scale in a given scheme for any given process. 
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BcGc(GeV- 2 ) 


No 
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0.63 


5.5 


0.417 


5.95-10" 4 


0.159 


0.46 


7.5 


0.558 


1.84-10" 6 


0.119 


0.43 


7.5 


0.565 


1.34-10 -6 


0.109 


0.54 


6.5 


0.484 


3.42-10" 5 


0.149 



TABLE II: Parameters of the b-CGC model obtained from fits to HERA data [lfj .The second row of parameters gives the best 
fit to HERA data. 
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FIG. 2: Saturation scales obtained from the IP-Sat(blue lines) and b-CGC models(red crosses) using first set of parameters in 
table. [J and the second set in table. [II] Left: impact parameter dependence of saturation scale different values of x. Right: x 
dependence of saturation scale for different values of impact parameter b. 



The NLO-BK Model 



The dipole cross-section is defined in eq. [2] as twice the forward scattering amplitude Af, which in the NLO-BK 
model [17J satisfies the equation, 



a/V(r±,y) 

dY 



= J dn K mo (r ± ,r 1 ,r 2 )x[Af(r 1 ,Y)+Af(r 2 ,Y)-Af(r ± ,Y)-Af(r 1 ,Y)Af(r 2> Y)] 



(15) 



where r 2 = rj_ — i"i and the NLO-kernel (where running coupling corrections are taken into account) is given by 



/CNLo(i\L,ri,r2) 
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1 f astf) 
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(16) 



This expression is based on considerable recent work to include running coupling corrections to the BK equation [311 
|32| . It should be noted that the expression does not include other next-to-leading log contributions to the kernel that 
have been computed recently [33| . Also of relevance to us is the assumption in the evolution equation in eq. [T5] that 
the dependence on the impact parameter and the dipole size factorize in the dipole amplitude as 



Af(r ± ,x,b ± ) = 2T{h ± )Af(r ± ,x). 



(17) 



The factorization here of impact parameter dependence and dipole size is very problematic conceptually and is an 
important limitation in applying these approaches to comparisons with data, except perhaps for final states that have 
limited sensitivity to the impact parameter dependence. How to include the impact pararmeter dependence in the 
BK / JIMWLK equations is an open question of great interest-we refer the reader to ref. 34| and references therein. 

The NLO-BK model was applied in refs. [35[ to a phenomenological study of the HERA data on the proton structure 
function F 2 . Here the impact parameter dependence is taken to be a step function. Two different models for the 
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initial condition were used in that work. The first is the GBW model 

N{r,Y = 0) = 1 - exp 

and the other is the MV model [2] 



N(r, Y = 0) = 1 - exp 



In 



rA 



QCD 



(18) 



(19) 



The initial condition for protons was determined from a global fit of F2 data in the work of [35|]. The fit parameters 
are summarized in the table HID The parameter C 2 here is a scale that controls the running of the coupling in 



i.e. 


o-p (fm 2 ) 


Q\ v (GeV 2 ) 


C 2 


7 


GBW 


3.159 


0.24 


5.3 


NA 


MV 


3.277 


0.15 


6.5 


1.13 



TABLE III: Parameters for the initial condition of the proton dipole cross section obtained in [351 ]. 

ref. [35[ . This parametrization was extended to nuclei and fit to fixed target e+A DIS data [24| ■ It was then applied 
to study long range correlations in A+A collisions [lij], single inclusive [HJ and double inclusive [22j distributions in 
deuteron-gold collisions at RHIC. In this work, only the MV parametrization will be considered. 



III. SINGLE INCLUSIVE DISTRIBUTION IN P+P COLLISIONS 

We shall now discuss the computation of the single inclusive hadron distribution in proton-proton collisions at the 
LHC obtained from the different dipole parametrizations discussed in the previous section. These parametrizations, 
when inserted in eqs. |4] and [3J respectively, give the single inclusive gluon distribution. While it might appear that 
all the parameters are fixed by the fits to the HERA data, a comparison with the single inclusive hadron distribution 
quickly makes clear that one requires further assumptions and parameters. In making comparisons to the collider 
inclusive hadron data, we shall discuss the sensitivity of our results to variations in these parameters. The additional 
assumptions and parameters are as follows: 



The dipole cross-sections are fit to HERA data for x < 0.01. One therefore needs to make an assumption for 
4>(x, fcj_, b±) for larger x > Xp = 0.01 values that kinematic regions of the proton-proton data are sensitive to. 
We use the parametrization [ljj 

( 1 — A ^ ( \ ^° 

(f>(x,kj_,b±_) = I :j— — - j ^—J (j)(x ,k±,b±_), x > x Q . (20) 

This parametrization of the large x unintegrated gluon distribution is motivated by quark counting rules (36j 
with fixed /3 = 4 and the parameter Ao which ranges from 0-0.2 in fits to the data. 

• The single inclusive gluon distribution in cq. [3] has a logarithmic infrared divergence which can be regulated 
either by putting a cutoff on lower limit of p± or replacing p± by m± — \Jp\+ m 2 where the mass is a free 
parameter. The single inclusive p± distribution is sensitive to the choice of m, which is fixed to be the same for 
data at all energy ranges. 

• Eq. [3] corresponds to the rapidity distribution of inclusive gluons at a fixed impact parameter. However, what 
is measured is the pseudo-rapidity distribution; the rapidity can be expressed in terms of the pseudo-rapidity 
most generally as 



y/ m 2 + p^ch 2 !] + p±shr] 
y/m 2 + p\ch 2 T] — p±shrj 



(21) 



The single inclusive distribution with respect to the pseudo-rapidity therefore contains a Jacobian from the 
conversion of the expression with respect to the rapidity. We choose for economy of parameters the mass term 
in eq. HT]to be the same as the one that regulates the infrared divergence. 
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• The minimum-bias single inclusive gluon distribution is obtained from the expression 



dN g fd*b ± £±- y (b ± ) 
d 2 p±dy J d 2 h± 



(22) 



In the IP-Sat model, the proton profile does not change with energy. The transverse diffusion of the proton, often 
termed Gribov diffusion [371 ] , is not fully accounted for by the diffusion of the unintegrated single inclusive gluon 
distribution because this growth does not automatically ensure the proper growth of the inelastic cross-section. 
Towards this end, we parametrize the maximum limit of b-integration to take the form b max . = &o + C ln(s). 
Fitting the available data on average dN/d?7 as a function of energy we can extract bo and C. For first set of 
parameters in table [Q the IP-Sat model gives &o=5.17 GeV -1 and C=0.19 with a choice of the mass term=0.4 
GeV used in the Jacobian. This value of bo is close to 2 & rms . in the IP-Sat model, where & rms . is the root mean 
square gluonic radius of the proton. The quantity irb"^^ , which is the denominator of eq. [55] can be interpreted 
as being closely related to the inelastic cross section contributing to particle production. Fig. [3] shows the 
variation of ^b 2 nax in the IP-Sat model as a function of collision energy. These numbers are in the ballpark of 
estimates of the inelastic cross-sections at the LHC |38|-thev are however significantly higher than the values at 
lower energies. This is because the numbers for & max . are extracted from fits of eq. [2"2l to data; it is therefore also 
sensitive to the uncertainties in the numerator of eq. 1221 A possible interpretation is that these uncertainties 
are larger at lower energies thereby leading to an overestimate of the inelastic cross-section. 

An additional point with regard to fig. [3] is that one observes a ~ 15% difference of 7r6^ ax at the highest LHC 
energies and may thereby hope to constrain the parameter m. However, due to the uncertainty in distinguishing 
the genuine inelastic cross-section from the NSD cross-section, it is unlikely at present that the 15% difference 
can be definitive in that regard. A similar form of b max when used to fit average dN/d?7 for b-CGC model gives 
C ~ 0. This suggests that because the impact parameter dependence of the dipole cross-section in the b-CGC 
model is tied in with its x dependence (see eq. II 2|) , the non-trivial relation of the two in this model may well 
approximate the physics of Gribov diffusion. 
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FIG. 3: Variation of the inelastic cross section with collision energy as obtained in the IP-Sat model for two different values of 
the infrared cut-off. See text for a detailed discussion. 

• The single inclusive hadron distribution is obtained by convolving eq. [22] with the fragmentation function for 
gluons into charged hadrons 1 , 



dN h 



d 2 p±dy 



dz dN g i n , 

. z 2 d 2 q±dy 



(23) 



1 We thank A. Dumitru for pointing out an error in this expression in an earlier version of the paper. After correcting the error, the 
resulting expression gives better agreement with data. 
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where D g ^h(z, fi 2 ) is chosen to be 6.05 z~°- 714 (l — z) 2 - 92 , corresponding to the LO parameter set of 39]. 
The lower limit of the integral is determined from the kinematic requirement that < 1. Our choice of 
fragmentation function is different from the fragmentation prescription of ref. [l2[ where the p± of gluons is 
scaled by p±/(z), with (z) = 0.5. We also note that we do not have additional scales such as intrinsic in our 
comparison to data, in contrast to ref. [l2[ . 



A. Results for the single inclusive hadron distribution 

Now that we have specified all the assumptions and parameters that go into eq. [23] we are ready to present the 
results of our comparisons of the different saturation models with the data on rapidity and transverse momentum 
distributions for a wide range of collider energies up to the highest present LHC energy of \/s = 7 TeV and make 
projections for ^/s = 10, 14 TeV. 




FIG. 4: dN/dr] obtained from comparing IP-Sat and b-CGC models to data from UA5 [!}, ALICE 0] and CMS The 
solid green band corresponds to uncertainties from different parameters; the dashed band is due to the variation of the choice 
of mass term in the Jacobian relating y to rj. The two curves at the top in both panels correspond to projections in the two 
models for *Js = 14 and 10 TeV respectively. 

Fig. 2] shows the pseudo-rapidity distribution obtained by integrating eq. [3] over p± for different mass terms in the 
Jacobian and parameters given in the tables Q] and |TT] for the IP-Sat and b-CGC dipole models respectively. The 
uncertainties corresponding to the choices of parametrizations and the infrared mass scale are shown in bands. In 
the IP-Sat model, the normalization is performed to the yfs dependence of dN/drj at r\ = as shown in fig. [6j-this 
fixes the two parameters in & max . we discussed previously. In the b-CGC model, because there is one less parameter 
(the coefficient of In s in 6 max . is zero) , the normalization can be performed to the rapidity distribution at one fixed 
energy. In the figure shown, the normalization is performed for y/s — 900 GeV; choosing a lower energy corresponds 
to a < 10% uncertainty in the overall normalization. 

We didn't use a fragmentation function in computing the pseudo-rapidity distributions because the rapidity distri- 
bution is vastly dominated by contributions below p± = 1 GeV, where fragmentation functions are likely not reliable. 
We have varied the mass term in the Jacobian corresponding to eq. [21] in the range 0.2-0.4 GeV, corresponding to 
an infrared scale of order Aqcd- (In each case we chose this mass term to be equal to the one we use to regulate the 
infrared divergence of the unintegrated gluon distribution in eq. 2]) . The effect of the extrapolation parameter Ao in 
eq. [20] is significant only at lower energies and and higher values of r\. The agreement with data of both models is 
quite good with the IP-Sat model providing a somewhat better agreement at the highest energies. 

In fig. [5] the corresponding p± distributions, with the previously specified fragmentation prescription, is shown. 
The IP-Sat model shows a poor agreement with data for the lower energies at high p±, as does the b-CGC model 
with some of the HERA parameter sets. At higher p±, the results are sensitive to physics at x > 0.01, which is 
parametrized very simply in the models. On the reverse side of the coin, one should anticipate a better agreement in 
the same p± window at higher energies. Indeed, a systematically better agreement is seen in both models with the 
data at the higher energies for all parameter sets. 
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FIG. 5: <m/d 2 p±dr) in the IP-Sat and b-CGC models. The solid green band corresponds to uncertainties from different 
parameters and the dashed band is due variation of the mass term in the Jacobian relating ri to y. The p± distribution is 
averaged over the rj range of ±2.4. The experimental data points are from CMS [55{ |. STAR [5b | and ATLAS [HtJ 



Fig. [6] shows the average value of dN/drj calculated at i]—0. The left plot of fig. [6] shows a fit of dN/dr] using 
different functional forms for the saturation scales. We considered both Q| and Q s / a s(Qs)- In the CGC framework, 
one would expect the latter. However, if the running is not significant in the energy range of interest, the former is 
also possible. Indeed, a good fit to the CMS data was obtained [44] with a simple form of the saturation scale a la the 
Golec-Biernat model 0- In our case, the comparison is made within the framework of the IP-Sat and b-CGC models 
which give better fits to the HERA data and are sensitive to the impact parameter profile of the gluon distribution 
in the proton. For the comparison with the LHC data in the left panel of fig. [Bl we use the value of Qs at the 
median value of s±=2 GeV -1 . The dependence of dN/drj on the purely Q s functional form is not very good, while 
the Qs/as{Qs) form does much better for the IP-Sat model. For the running of as, we chose Q s (s±) at s±=0 to 
restrict its running to as below 0.5. Fig. [6] (right panel) shows by way of comparison, a comparison of dN/dr) at 
7] = as a function of ^/s to IP-Sat and b-CGC models. In the IP-Sat model, a good fit is ensured because a fit to 
this energy dependence is what determines the parameters of 6 ma x.-see fig. [3] and related discussion. In the b-CGC 
model, the curve is a prediction and is seen to be a very good fit to the data. 

The energy dependence of (p±) is shown in fig. [71 In the left panel, it is shown that in this case one obtains a good 
linear dependence of (p±) on Qs in both the IP-Sat and b-CGC models as the cm energy is varied. This is as seen 
previously [Hj|. The right plot shows (p±) versus ^/s computed in the IP-Sat and b-CGC models. Here one sees that 
the results are quite sensitive to choice of the infrared cut-off. 

Fig. [5]shows the r\ and p± distributions computed in the NLO-BK model. Only MV initial conditions are considered. 
In this model, the impact parameter dependence of the inclusive gluon multiplicity is given by 

Tib) = — ±— 8 (6 max . - b) . 

u max. 

Therefore dN g (b)/dr] = dN g /dr] Tib), where 6 max . is a parameter that can be absorbed in the normalization. If no 
dependence of 6 max . on y/s is assumed, the model considerably overestimates the single inclusive data at LHC energies. 
This is likely a consequence of the fact that the NLO-BK doesn't take into account the impact parameter dependence 
of the saturation scale. In this case, the values for the inclusive gluon multiplicity correspond to the values for the 
zero impact parameter, which is considerably higher than the minimum bias values. The agreement with data is 
improved considerably by allowing & max . to depend on y/s. The denominator Trti? nax therefore provides an energy 
dependent normalization. The same approach discussed previously for the IP-Sat model is employed to extract the 
yfs dependence of 6 max . chosen to be of the form 6 max . = bo + C ln(s). From a fit to the average dN/dr/ as a function 
of energy, one obtains 6q ~ 5.6 — 7.55 GeV -1 and C ~ 0.23 — 0.46 depending on the choice of infrared cut-off. 
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FIG. 6: Average dN/dq\ in the IP-Sat and b-CGC models. Left: Data plotted as a function of saturation scales (at median 
impact parameter 6 = 2 GeV -1 ) for both the models determined from the HERA data. Right: Average dN/dr) at r\ = from 
the fcj_-factorized expression in eq. [23] from IP-Sat (solid green) and b-CGC (dashed) models. See text for further discussion. 
Experimental data points are from Ref.[H|], & 03, HI], HI- 
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FIG. 7: Average p± obtained from IP-Sat and b-CGC models compared to data. Left: Function of saturation scales fitted 
for both the models. Right: Average p±_ from the fcx-factorized expression in eq. [23] The thin (colored) bands correspond to 
uncertainties arising from different parameters in tableland tableland mass term=0.2 GeV. The thick (gray) band shows the 
sensitivity to variation of mass term in the range 0.2-0.3 GeV (with lower mass corresponding to lower (p±) ). Experimental 
data points are from Ref. [H [13 & 



IV. MULTIPLICITY DISTRIBUTION 



There are several sources of multiplicity fluctuations in high energy hadronic collisions. These can arise from 
fluctuations in the number of wee partons, in their distribution with impact parameter and their distribution in 
rapidity [45| . In this paper, we will consider particle production in a relatively small rapidity window (parametrically 
of order Arj < 1/as), so fluctuations in rapidity will not be an important source of fluctuations. Let us first consider 
fluctuations in multiplicity for a fixed impact parameter. In this case, the CGC framework allows for a systematic 
treatment of inclusive multi-particle production [47j in the Glasma |4g|. The largest contribution to multi-particle 
production comes from diagrams that appear superficially disconnected, but are connected by averaging over color 
correlations in an event and over all events. This is the formal basis of the Glasma flux tube picture [481 1 . an d was 
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FIG. 8: Pseudo-rapidity and p± distribution in the NLO-BK model compared to data. The uppermost two plots in the left 
panel correspond to predictions for yfs = 14, 10 TeV with m = 0.4 GeV. The p± distribution is averaged over the n range of 
±2.4. The band corresponds to the variation m = 0.2-0.4 GeV in the Jacobian relating r\ to y. 



previously used to compute backward- forward correlations [43, [50] , two particle correlations [24|, [48[ , three particle 
correlations [5l|, and n- particle correlations (52[. It is this last computation that will concern us here. The n- particle 
correlations obtained by averaging over color sources in the CGC picture are those that would be generated by the 
negative binomial distribution [52j j . To the best of our knowledge, this derivation of the negative binomial distribution 
as arising from particle production from Glasma flux tubes is the first such ab initio derivation in a QCD framework 2 
The negative binomial distribution is given by 



pNB f - k) _ r ( fc + n ) ^ , ■>,, 

" 1 ' >~ r(fc)i> + i) (n + k)^ ■ 



The distribution is characterized by two parameters, the mean multiplicity n and the parameter k. As is well known, 
in the limit k — > oo, this distribution reduces to the Poisson distribution. In the limit k — > 1, one obtains the Bose- 
Einstcin distribution. The variance of the distribution is given by a 1 = n 2 — n 2 = n + n 2 /k. In the Glasma flux tube 
approach, k is not an arbitrary parameter; instead, it is computed to be 

where £ is a dimensionless non-perturbative parameter which will be discussed shortly and S± is the overlap area of 
the two hadrons. 

As we mentioned previously, the negative binomial distribution is obtained at a fixed impact parameter, so n = n(b) 
and k = k(b). In particular, the latter parameter must be interpreted as being proportional to the number of flux 
tubes (or interacting "hot spots") S_\_/1/Q 2 S at a given impact parameter. Because Q 2 S grows with energy, k has a 
very particular energy dependence. Here, the parameter k in eq. [35] is determined as follows. For a given impact 
parameter, we define 



Q%S X = J d 2 x±Q 2 { Xl _), 



where the integral on the r.h.s is performed over the overlap area of the two protons at a given parameter. At 

)^,Qf}, where Qj, and Qf 



each given transverse position x± in the overlap area, we choose Qs(x±) = min. {Qg, Qg }, where Qg and Q$ are 



2 The negative binomial distribution has of course been known for a long time [641 to provide good fits to collider p+p data. More recently 
shown to describe the multiplicity distributions in A+A collisions at RHIC [65j |. However, a microscopic derivation of this distribution 
from QCD was previously lacking. 
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respectively the saturation scales of the two colliding protons at that x±. This choice is motivated by the fact that 
the inclusive multiplicity of produced gluons is much more sensitive to the smaller of the two saturation scales [53| . 
In fig. [9] (left) we plot the saturation scale as a function of y/s for different impact parameters. The right plot has on 
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FIG. 9: Left: Variation of the saturation scale with the collision energy. Right: Variation of Q%S± for different impact 
parameters and the cm. energy of the collision. The solid lines are for the IP-Sat model and crosses for the b-CGC model, in 
each case for the parameters providing the best fit to the HERA data. 



the y-axis the quantity that appears in eq. [55] determined by the procedure we described. We observe that a stronger 
dependence of this quantity is seen for the b-CGC model. With the exception of the parameter £, we have everything 
necessary to compute eq. [24] at a given impact parameter. 

Fluctuations in impact parameter are treated as follows. The overlap function for two protons at a given impact 
parameter (see fig. Q] (right)) can be expressed as 



T PP {h±_) = J dVT p (sx)T p (s ± -bj 



(26) 



where T p is given for instance in the IP-Sat model by eq. [TT] Our knowledge of the HERA diffractive data therefore 
allows one to compute T pp in the saturation models. The probability distribution for an inelastic collision at a given 
impact parameter is given in impact parameter eikonal models as [3a . [ill . |42[ 



jpcik. 
ar incl. 



1 - CXp (-CTggTpp) 



d 2 b ± J d 2 b ± (1 - exp (-OggTpp)) 



(27) 



In general, a m is an energy dependent quantity estimated to be the elementary cross-section for gluon-gluon scattering. 
Alternately, in our framework, an estimate for this quantity in our framework involving no additional parameters is 



dP 



dip. 



dN g 
dy 



(bj 



d 2 b. 



(28) 



This expression of course gives unity when both sides are integrated over impact parameter. The result for 
2tt6 dP^ h /d 2 b± at a fixed = 900 GeV is shown in fig [10] and is a sharply peaked distribution at ~ 3 GeV . 
The distribution plotted is insensitive to yfs and to the infrared scale to. We should note here that the inclusive 
multiplicity at a given impact parameter is given by fi(h±) = Pi + 2 Pi + 3 P3 + • • • , while the inelastic probability 
is given by Pinci. = Pi + P% + P3 + ■ ■ ■ . In general, the two are of course not the same, so using eq. [55] to estimate 
the inelastic probability is strictly not correct. Our justification here is that eq. 1281 has the right qualitative behavior 
and involves one less parameter than eq. [27] It is very desirable to obtain a better estimate for eq. [28] For a detailed 
discussion of the problem of computing multi-particle production in the presence of strong time dependent sources, 
see ref. EH- 

With the stated assumptions (and caveats) we are now in a position to compute the probability distribution as 
a function of multiplicity. By convolving the probability distribution for producing n particles at a given impact 
parameter (eq. [24]) with the probability for an inelastic collision at that impact parameter (from eq. l28|). one obtains 
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FIG. 10: Probability distribution as a function of impact parameter for an inelastic collision computed using eq. [28] for the 
b-CGC and IP-Sat models. 



The results for this quantity are shown in fig. 1111 Note that since the input here is the average inclusive multiplicity 
at a given impact parameter (as opposed to the minimum bias inclusive multiplicity that was compared to data), 
this quantity needs to be normalized as well. We do so by fitting the multiplicity distribution corresponding to the 
lowest energy UA5 data set at y/s = 200 GeV for the normalization of n. In the b-CGC model, a fit of the overall 
normalization to the single inclusive minimum bias distribution at a given energy gives the same value ( to an accuracy 
of < 5 %) as that obtained for the same quantity if fit to the multiplicity distribution instead. In the IP-Sat model, one 
obtains the same normalization constant if 6 max . = 2 b Tms .. If one chooses the form 6 max . = &o + Cm(s) we described 
previously, there is a 60% discrepancy between the two choices of fixing the normalization. The agreement between 
data and model is remarkably good for IP-Sat distribution for all energies and rapidity cuts 3 . For the b-CGC model, 
the agreement with data for |?7| < 0.5 is quite good for the 2.36 TeV ALICE data [6(| but shows deviations for other 
energies at the highest multiplicities 4 . 

An important point regarding the comparison of the models to data in fig. [11] concerns the magnitude of the 
parameter £ in eq. [55] which is fit to the data. In the b-CGC model, it is extracted to be 0.25 and it is 0.35 in the 
IP-Sat model 5 . This quantity was recently computed non-perturbatively ab initio by solving the Yang-Mills equations 
numerically (6(| for two gluon correlations from gauge fields generated in the collision of two dense color sources (6(| . 
The results of the numerical computation vary depending on parameter choices in the range C ~ 0.3-1.5-the essential 
point is that the result is a number of order unity. It is very encouraging that the values extracted from the data 
for ( in the saturation models are numbers of this order 6 , and suggests that careful analysis of the data can help to 
extract non-trivial information about the screening of Glasma flux tubes in p+p collisions. 



3 The results are insensitive to the choice of the infrared cut-off m. 

4 Our computation of the saturation scale is performed for r\ = while the data we compare with are averaged over < 0.5 and \f)\ < 1. 
This discrepancy may lead to small corrections to our results. 

5 We have not attempted any fine tuning of our fits at this stage. A good agreement can be obtained for the 2.36 TeV and 7 TeV data in 
the b-CGC model by decreasing f by 20% but this gives a significantly worse fit at lower energies. However, because the model's validity 
is questionable at lower energies, and because of other sources of uncertainty we have articulated, the b-CGC model still provides a 
"competitive" mechanism for multi-particle production. 

6 Nothing in principle would have prevented the number for f extracted from the data from being orders of magnitude off, in which case 
the present analysis would be invalidated completely. 



the expression 




(29) 
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FIG. 11: Probability distribution of the gluon multiplicity computed in the IP-Sat model (parameter set I) and the b-CGC 
model (parameter set II) compared to the UA5[H^] and ALICE [6C| data for different r\ ranges. Left: Multiplicity distribution 
for the r\ range ±0.5. Right: Multiplicity distribution for the r\ range ±1. 



V. SUMMARY AND OUTLOOK 



In this work, we extracted impact parameter dependent unintegrated gluon distributions from fits to the HERA 
inclusive and exclusive data in the IP-Sat, b-CGC and NLO-BK models. These models all implement the physics of 
saturation but differ in their dynamical assumptions. The impact parameter distributions, in the k± factorization 
formalism allow one to compute single inclusive rapidity and p± distributions in p+p collisions at the LHC. These 
give quite reasonable agreement with the LHC data. These impact parameter dependent distributions also allowed us 
to compute the multiplicity distribution, which in the Color Glass Condensate/ Glasma formalism is predicted to be a 
negative binomial distribution with particular values for the parameter controlling the width of the distribution. We 
observed that the multiplicity distributions are well reproduced in this framework with numbers for the parameter 
consistent with values extracted from numerical computations of Yang-Mills equations for the collision of dense 
color sources. These results suggest that particle emission from Glasma flux tubes generated in collisions of "hot 
spots" of size 1/Qs are a strong candidate for generating the multi-parton correlations underlying the multiplicity 
distribution. Because these Glasma flux tubes generate long range rapidity correlations, our results are also consistent 
with a recent claim I68ll that Glasma flux tubes generate the near side ridge seen in high multiplicity events by the 
CMS collaboration [67] . Because our analysis allows a more careful treatment of the contribution of various impact 
parameters to the multiplicity distribution, it will allow more quantitative and systematic comparisons to phenomena 
seen only in high multiplicity cuts in the LHC data. 

There are several caveats in this analysis that must be noted. Firstly, the k± factorization formalism for inclusive 
multiplicity distributions is rather fragile for k± < Q$ because it does not include multi-parton rescatterings [69l [7fj|. 
Significant improvements of the k± factorization formalism are feasible in the CGC framework if presently cumbersome 
to implement in phenomenological analyses. Similarly, NLO dipole computations are becoming available albeit a 
consistent NLO dipole analysis still remains to be developed. Both these developments suggest that a combined 
quantitative study of the collective QCD dynamics of saturation in DIS and hadronic collisions is feasible in future 
along the lines discussed in this work. 



Acknowledgements 

R.V was supported by the US Department of Energy under DOE Contract No.DE-AC02-98CH10886. We thank 
Kevin Dusling and Frangois Gelis for a careful reading of the manuscript and especially Adrian Dumitru for pointing 
out an error in a previous version of the manuscript. We gratefully acknowledge useful conversations with Javier 



15 



Albacete, Guillaume Beuf, Subhasis Chattopadhyay, Tuomas Lappi, Larry McLerran, Zhangbo Kang and Feng Yuan. 



[1] L.V. Gribov, E.M. Levin, M.G. Ryskin, Phys. Rept. 100, 1 (1983); A.H. Mueller, J-W. Qiu, Nucl. Phys. B 268, 427 
(1986). 

L.D. McLerran, R. Venugopalan, Phys. Rev. D 49, 2233 (1994); ibid. 49, 3352 (1994); ibid. 50, 2225 (1994). 
Jalilian-Marian, A. Kovner, A. Leonidov, H. Weigert, Nucl. Phys. B 504, 415 (1997); ibid., Phys. Rev. D 59, 014014 
(1999); E. Iancu, A. Leonidov, L.D. McLerran, Nucl. Phys. A 692, 583 (2001); E. Ferreiro, E. Iancu, A. Leonidov, L.D. 
McLerran, Nucl. Phys. A 703, 489 (2002). 



E. Iancu, R. Venugopalan, Quark Gluon Plasma 3, Eds. R.C. Hwa, X.N. Wang, World Scientific, hep-ph/0303204 

H. Weigert, Prog. Part. Nucl. Phys. 55, 461 (2005); F. Gelis, E. Iancu, J. Jalilian-Marian, R. Venugopalan. larXiv:1002.0333 
L. D. McLerran, R. Venugopalan, Phys. Rev. D 59, 094002 (1999); R. Venugopalan, Acta Phys. Polon. B 30, 3731 (1999) 
N. N. Nikolaev, B. G. Zakharov, Z. Phys. C 49, 607 (1991); A. H. Mueller, Nucl. Phys. B 335, 115 (1990). 

K. J. Golec-Biernat, M. Wusthoff, Phys. Rev. D 59, 014017 (1998); K. J. Golec-Biernat, M. Wusthoff, Phys. Rev. D 60, 
114023 (1999) 

I. Balitsky, Nucl. Phys. B 463, 99 (1996); Yu.V. Kovchegov, Phys. Rev. D 61, 074018 (2000). 
J. P. Blaizot, F. Gelis, R. Venugopalan, Nucl. Phys. A 743,57 (2004). 

MA. Braun, Phys. Lett. B 483, 105 (2000). 

F. Gelis, A. M. Stasto, R. Venugopalan, Eur. Phys. J. C 48, 489 (2006). 
E. Levin, A. H. Rezaeian, Phys. Rev. D 82, 014022 (2010). 

E. Iancu, K. Itakura, S. Munier, Phys. Lett. B 590, 199 (2004). 
H. Kowalski, L. Motyka, G. Watt, Phys. Rev. D 74, 074016 (2006) 

G. Watt, H. Kowalski, Phys. Rev. D 78, 014016 (2008) 

H. Kowalski, D. Teaney, Phys. Rev. D 68, 114005 (2003). 

J. L. Albacete, Y. V. Kovchegov, Phys. Rev. D 75, 125021 (2007). 

J. R. Forshaw, R. Sandapen, G. Shaw, JHEP 0611, 025 (2006); B. Z. Kopeliovich, J. Raufeisen, A. V. Tarasov, Phys. Rev. 
C 62, 035204 (2000); . Gotsman, E. Levin, M. Lublinsky, U. Maor, Eur. Phys. J. C 27, 411 (2003); G. Soyez, Phys. Lett. 
B 655, 32 (2007). 

D. Kharzeev, Yu. Kovchegov, K. Tuchin, Phys. Rev. D 68, 094013 (2003); A. Dumitru, A. Hayashigaki, J. Jalilian-Marian, 
Nucl. Phys. A 770, 57 (2006). 

V. P. Goncalves, M. S. Kugeratski, M. V. T. Machado, F. S. Navarra, Phys. Lett. B 643, 273 (2006). 

J. L. Albacete, C. Marquet, Phys. Lett. B 687, 174 (2010). 

J. L. Albacete, C. Marquet, Phys. Rev. Lett. 105, 162301 (2010). 

H. Kowalski, T. Lappi, R. Venugopalan, Phys. Rev. Lett. 100, 022303 (2008). 
K. Dusling, F. Gelis, T. Lappi, R. Venugopalan, Nucl. Phys. A 836, 159 (2010). 
T. Rogers, V. Guzey, M. Strikman, X. Zu, Phys. Rev. D 69, 074011 (2004). 

J. Bartels, K. J. Golec-Biernat, H. Kowalski, Phys. Rev. D 66, 014001 (2002). 
S. Chekanov et al. [ZEUS Collaboration], Eur. Phys. J. C 21, 443 (2001). 

C. Adloff et al. [HI Collaboration], Eur. Phys. J. C 21, 33 (2001). 

S. Chekanov et al. (ZEUS Collaboration), Eur. Phys. J. C 24, 345 (2002), Nucl. Phys.B695, 3 (2004). 

A. Aktas et al. (HI Collaboration), Eur. Phys. J. C 46, 585 (2006). 

I. Balitsky, Phys. Rev. D 75, 014001 (2007). 

Y. V. Kovchegov, H. Weigert, Nucl. Phys. A 784, 188 (2007). 

I. Balitsky, G. A. Ch irilli, Phys. Rev. D 77, 014019 (2008); G. Beuf. larXiv: 1008.0498 1 [hep-ph] . 
J. Berger, A. Stasto. larXiv:1010.067T1 [hep-ph], 

J. L. Albacete, N. Armesto, J. G. Milhano , C. A. Salgado, Phys. Rev. D 80, 034031 (2009). 

V. A. Matveev, R. M. Muradian, A. N. Tavkhelidze, Lett. Nuovo Cim. 7, 719 (1973); S. J. Brodsky, G. R. Farrar, Phys. 

Rev. Lett. 31 1153 (1973). 

V. N. Gribov, |arXiv:hep-ph/0006158| 

D. d'Enterria, G. K. Eyyubova, V. L. Korotkikh, I. P. Lokhtin, S. V. Petrushanko, L. I. Sarycheva, A. M. Snigirev, Eur. 
Phys. J. C 66, 173 (2010). 

B. A. Kniehl, G. Kramer, B. Potter, Nucl. Phys. B 582, 514 (2000). 

S. Eidelman et al. [Particle Data Group Collaboration], "Review of particle physics", Phys. Lett. B592 1 (2004). 

L. Frankfurt, M. Strikman, C. Weiss, Ann. Rev. Nucl. Part. Sci. 55, 403 (2005). 

M. L. Miller, K. Reygers, S. J. Sanders, P. Steinberg, Ann. Rev. Nucl. Part. Sci. 57, 205 (2007). 

F. Gelis, R. Venugopalan, Nucl. Phys. A 776, 135 (2006); ibid. A 779, 177 (2006). 
L. McLerran, M. Praszalowicz, Acta Phys. Polon. B 41, 1917 (2010). 

H. I. Miettinen, J. Pumplin, Phys. Rev. D 18, 1696 (1978). 
T. Lappi, L.D. McLerran, Nucl. Phys. A 772, 200 (2006). 

F. Gelis, T. Lappi, R. Venugopalan, Phys. Rev. D 78, 054019 (2008); ibid., 78, 054020 (2008); ibid., 79, 094017 (2009). 
A. Dumitru, F. Gelis, L. McLerran, R. Venugopalan, Nucl. Phys. A 810, 91 (2008). 



[49] N. Armesto, L. McLerran, C. Pajares, Nucl. Phys. A 781, 201 (2007). 

[50] T. Lappi, L. McLerran, Nucl. Phys. A 832, 330 (2010). 

[51] K. Dusling, D. Fernandez-Fraile and R. Venugopalan, Nucl. Phys. A 828, 161 (2009) [arXiv:0902.4435l [nucl-th]] 
[52] F. Gelis, T. Lappi, L. McLerran, Nucl. Phys. A828 (2009) 149. 

[53] A. Dumitru, L. D. McLerran, Nucl. Phys. A 700, 492 (2002). 

[54] ALICE Collaboration, Eur. Phys. J. C 65 (2010) 111 |arXiv:0911.5430] ; larXiv: 1004.30341 
[55] CMS Collaboration, Phys Rev Lett 105,022002(2010)/ 

[56] STAR Collaboration, Phys. Rev. Lett. 91, 172302 (200 3), Phys. Rev . C 79, 034909 (2009). 
[57] ATLAS Collaboration, Phys. Lett. B688 (2010) 21 [arXiv: 1003.3124] . 
[58] UA5 collaboration, Z. Phys. C33 (1986) 1. 

[59] UA5 Collaboration, Z. Phys. C43 (1989) 357. 

[60] ALICE Collaboration, Eur.Phys.J.C68:345-354,2010 |arXiv: 1004.3514] . 

[61] UA1 collaboration, A. M. Rossi et al., Nucl. Phys. B84 (1975) 269. 

[62] CDF collaboration, Phys. Rev. D41 (1990) 2330; Phys. Rev. Lett. 61 (1988) 1819. 

[63] E735 Collaboration, T. Alexopoulos et al., Phys. Rev. Lett. 60, 1622 (1988). 

[64] E. A. De Wolf, I. M. Dremin, W. Kittel, Phys. Rept. 270, 1 (1996). 

[65] PHENIX collaboration, A. Adare et. al., Phys. Rev. C78 (2008) 044902. 

[66] T. Lappi, S. Srednyak, R. Venugopalan, JHEP 1001 066 (2010). 

[67] CMS Collaboration, I arXiv: 1009.41221 

[68] A. Dumitru, K. Dusling, F. Gelis, J. Jalilian-Marian, T. Lappi, R. Venugopalan, arXiv: 1009.5295 [hep-ph]. 
[69] A. Krasnitz, R. Venugopalan, Nucl. Phys. B 557, 237 (1999); A. Krasnitz, Y. Nara, R. Venugopalan, Phys. Rev. Lett. 
192302 (2001). 

[70] T. Lappi, Phys. Rev. C 67, 054903 (2003); T. Lappi, Eur. Phys. J. C 55, 285 (2008). 



